clear
set more off
capture postclose bskeep

local h_or_o_2012 = 0

use opioids_by_month, clear

drop if year<2004

* construct the dummies used in analysis
gen trend=12*(year-2004)+month
gen heroindr=(heroin_death/(population/100000))
gen opioidr=(opioid_deaths/(population/100000))
gen h_o_deathr=(h_o_deaths/(population/100000))
global mint=12*(2009-2004)+1
global maxt=12*(2012-2004)+12

preserve
use ../factors/pred_deaths.dta, clear
collapse (sum) h_hat_spec, by(year month)
sort year month
tempfile measured
save `measured'
	
restore

sort year month
merge 1:1 year month using `measured'
drop _merge
gen hhatdr = h_hat_spec/(population/100000)


tempfile main bootsave 


if `h_or_o_2012' == 1 {
	qui save `main' , replace 
	postfile bskeep trend sseh teststat pvalp r2 using sense_horo, replace

	
	global maxt=12*(2011-2004)+12
	forvalues b = $mint/$maxt { 
		display `b'
		use `main', clear 
		gen z=trend-`b'
		gen post=trend>=`b'
		gen zpost1=post*z
		gen zpost2=zpost1*zpost1
		gen zpre1=(1-post)*z
		gen zpre2=zpre1*zpre1
		
		** subset of years used in heterogeneity analysis
		reg h_o_deathr zpre1 zpre2 zpost1 zpost2 if year>=2004 & year<=2012
		local sseh=e(rss)
		test (zpre1=zpost1) (zpre2=zpost2)
		local teststat = r(F)
		local pvalp = r(p)
		local spec = 1
		post bskeep (`b') (`sseh') (`teststat') (`pvalp') (e(r2))
	}
	postclose bskeep
	clear
	use sense_horo, clear
	sort sseh
	gen keeper = _n==1
	list
	sort trend
	tempfile bests 
	save `bests'

	use `main', clear
	sort trend
	merge 1:1 trend using `bests'

	egen maxer = max(teststat)
	egen breaker = max( (teststat==maxer)*trend )
	gen post = trend>=breaker
	gen zpost1 = post*(trend-breaker)
	gen zpre1 = (1-post)*(trend-breaker)
	gen zpost2 = zpost1^2
	gen zpre2 = zpre1^2
	reg h_o_deathr zpre1 zpre2 zpost1 zpost2 if year>=2004 & year<=2012
	test (zpre1=zpost1) (zpre2=zpost2)
	predict yhatter
	
	keep year month trend h_o_deathr teststat yhatter
	keep if year>=2004 & year<=2012
	export excel using ./h_or_o_drate_spline.xls, replace
	
	*list year month trend teststat

}


